This is a step-by-step guide on how to perform Remote Sensing Image Analysis and Change Detection (NDVI) in R.

Tools:

Background:

Normalized Difference Vegetation Index (NDVI) quantifies vegetation by measuring the difference between near-infrared (which vegetation strongly reflects) and red light (which vegetation absorbs). NDVI always ranges from -1 to +1.

For example, when you have negative values, it’s highly likely that it’s water. On the other hand, if you have an NDVI value close to +1, there’s a high possibility that it’s dense green leaves.

But when NDVI is close to zero, there are likely no green leaves and it could even be an urbanized area.

NDVI is the most common index that analysts use in remote sensing. (Source: gisgeography.com)

Install packages/Load Package 

library(raster)
## Loading required package: sp
library(sf)
## Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE
library(sp)
library(mapview)

Set Working Directory

It is always best to set working directories as all files will be written or obtain from this location.

wd <- "C:\\path\\to\\your\\directory\\L2A_T16QCE_A026335_20220322T162650_2022-03-22_con"
setwd(wd)

Load Raster Bands to RStudio

Spectral remote sensing data are collected by powerful camera-like instruments known as imaging spectrometers. Imaging spectrometers collect reflected light energy in “bands.”

A band represents a segment of the electromagnetic spectrum. You can think of it as a bin of one “type” of light.

Atmospheric transmission of the different bands from Earth.

#The file name is constant so we can copy the string and paste the different band numbers 
file_name = "\\RT_T16QCE_A026335_20220322T162650_B"
# blue band
b2<- raster(paste0(wd,file_name,"02.tif", sep=""))
# green band
b3<- raster(paste0(wd,file_name,"03.tif", sep=""))
# red band
b4<- raster(paste0(wd,file_name,"04.tif", sep=""))
#veg red Edge = the region where the spectral reflectance of green vegetation rises rapidly within a certain band range
b5<- raster(paste0(wd,file_name,"05.tif", sep=""))
b6<- raster(paste0(wd,file_name,"06.tif", sep=""))
b7<- raster(paste0(wd,file_name,"07.tif", sep=""))
# NIR = near infrared
b8<- raster(paste0(wd,file_name,"08.tif", sep=""))
# vegetion red edge
b8A<- raster(paste0(wd,file_name,"8A.tif", sep=""))
#short wave infrared
b11<- raster(paste0(wd,file_name,"11.tif", sep=""))
b12<- raster(paste0(wd,file_name,"12.tif", sep=""))

Stack the raster bands and explore the data.

bzRGB <- stack(list(b4, b3, b2))
bzRGB
## class      : RasterStack 
## dimensions : 10980, 10980, 120560400, 3  (nrow, ncol, ncell, nlayers)
## resolution : 10, 10  (x, y)
## extent     : 3e+05, 409800, 1890240, 2000040  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs 
## names      : RT_T16QCE_A026335_20220322T162650_B04, RT_T16QCE_A026335_20220322T162650_B03, RT_T16QCE_A026335_20220322T162650_B02

View the histogram of the red, green and blue bands

par(mfrow=c(2,2))
hist(b4, col='#FF0000', main = "Band 04: Red")
hist(b3, col ='#50C878', main = "Band 03: Green") 
hist(b2, col='#0000FF', main = "Band 02: Blue")

View the image in human readable true color composite.

plotRGB(bzRGB, axes = TRUE, stretch = "lin", main = "Sentinel RGB True Color Composite")

Calculate the NDVI

The Normalized Difference Vegetation Index takes into account the amount of near-infrared (NIR) reflected by plants. It is calculated by dividing the difference between the reflectances (Rho) in the near-infrared and red by the sum of the two. NDVI values typically range between negative one (surface water) and one (full, vibrant canopy). Low values (0.1 – 0.4) indicate sparse canopies, while higher values (0.7 – 0.9) suggest full, active canopies.

#Calculate the Normalized Difference Vegetation Index
ndvi <- (b8 - b4) / (b8 + b4)

#View the NDVI
plot(ndvi, rev(terrain.colors(10)), main = 'Sentinel-2 NDVI - Belize Coast')

See the distribution of the NDVI Values.

hist(ndvi,
     main = "Distribution of NDVI values",
     xlab = "NDVI",
     ylab= "Frequency",
     col = "aquamarine3", #color in R 
     xlim = c(-0.5, 1),
     breaks = 30,
     xaxt = 'n') # stat frequency 
axis(side = 1, at = seq(-0.5,1, 0.05), labels = seq(-0.5,1, 0.05))

We are reclassifying our object and making all values between negative infinity and 0.4 be NAs

veg <- reclassify(ndvi, cbind(-Inf, 0.4, NA))

We plot the Vegetation Cover.

plot(veg, main = 'Vegitation Cover')

We can also create an interactive map to see the data.

mapviewOptions(basemaps.color.shuffle = FALSE) # change default enhance contrast based on layer color
mapviewOptions(basemaps =c("Esri.WorldImagery", "CartoDB.DarkMatter", "CartoDB.Positron", "OpenStreetMap", "OpenTopoMap")) #change order of basemap
map<- mapview(veg) #create a mapview object
map2<- mapview::viewRGB(bzRGB, r = 1, g = 2, b = 3) # true-color of raster stack
map3 <- map2 + map #add maps together
map3 # view the map


We can also write raster to file and open in other Desktop GIS application like QGIS.

#writeRaster(x = veg,
            # where your file will go - update with your file path!
#            filename="C:\\path\\to\\your\\directory\\veg_2022.tif",   
#            format = "GTiff",                  # save as a tif
#            datatype = 'INT2S')
–END–